#include <stdio.h>
#include <gmp.h>
#include <math.h>
#include <pthread.h>

#define IT 10000
#define PREC 1000

// funções para thread
void *calc_a(void *t);
void *calc_b(void *t);
void *calc_t(void *t);
void *calc_p(void *t);

// variaveis globais
mpf_t a_prev, a_next;
mpf_t b_prev, b_next;
mpf_t t_prev, t_next;
mpf_t p_prev, p_next;


// iniciando o programa
int main(int argc, char *argv[]) {

	// verefica se a entrada esta correta
	if (argc != 2) {
		printf("usage: ./gauss_legendre (0 = sequential / 1 = parallel)\n");
		return -1;
	}
	
	// seta PREC como precisão
	mpf_set_default_prec(PREC);

	// inicializa as variaveis gmp (big num)
	mpf_init(a_prev); mpf_init(a_next);
	mpf_init(b_prev); mpf_init(b_next);
	mpf_init(t_prev); mpf_init(t_next);
	mpf_init(p_prev); mpf_init(p_next);

	// cria variavel pi do tipo gmp
	mpf_t pi;

	// inicializa pi
	mpf_init(pi);

	int i;
	double d = 1;


	// INICIO DO CALCULO
	// Aqui a primeira parte do calculo como não é nenhuma interação,
	// apenas algumas linhas a serem executadas (que são necessárias
	// para começar a interação), foi feito em sequencial para serem usadas
	// depois tanto com o processo em paralelo como em sequencial

	// a_next = d
	mpf_set_d(a_next, d);
	// p_next = d
	mpf_set_d(p_next, d);

	// 1/raiz(2)
	d = 1 / sqrt(2);
	
	// b_next = d
	mpf_set_d(b_next, d);

	d = 0.25;
	// t_next = d
	mpf_set_d(t_next, d);


	// INICIO DO PROCESSO EM SEQUENCIAL
	if (argv[1][0] == '0') {

		// começa a interação de calculos em sequencial
		for(i=0; i<IT; ++i) {
			mpf_set(a_prev, a_next);			
			mpf_set(b_prev, b_next);			
			mpf_set(t_prev, t_next);			
			mpf_set(p_prev, p_next);			

			// a_next = (a_prev + b_prev) / 2;
			mpf_add(a_next, a_prev, b_prev);
			mpf_div_ui(a_next, a_next, 2);

			// b_next = (a_prev * b_prev) ^ 1/2
			mpf_mul(b_next, a_prev, b_prev);
			mpf_sqrt(b_next, b_next);

			//t_next = t_prev - p_prev (a_prev - a_next) ^ 2
			mpf_sub(t_next, a_prev, a_next);
			mpf_pow_ui(t_next, t_next, 2);
			mpf_mul(t_next, t_next, p_prev);
			mpf_sub(t_next, t_prev, t_next);

			// p_next = 2 * p_prev
			mpf_mul_ui(p_next, p_prev, 2);

		}

	// INICIO DO PROCESSO EM PARALELO
	} else {

		// cria threads
		pthread_t a_next_t;
		pthread_t b_next_t;
		pthread_t t_next_t;
		pthread_t p_next_t;

		// começa a interação em paralelo
		for(i=0; i<IT; ++i) {
			mpf_set(a_prev, a_next);			
			mpf_set(b_prev, b_next);			
			mpf_set(t_prev, t_next);			
			mpf_set(p_prev, p_next);			

			// thread para calcular 'a', 'b' e 'p'
			pthread_create( &a_next_t, NULL, calc_a, NULL);
			pthread_create( &b_next_t, NULL, calc_b, NULL);
			pthread_create( &p_next_t, NULL, calc_p, NULL);

			// espera thread de calculo de 'a' acabar
			// pois a próxima instrução depende de 'a_next'
			pthread_join(a_next_t, NULL);
			
			// calcula t_next
			pthread_create( &t_next_t, NULL, calc_t, NULL);

			// espera todas as threads finalizarem
			// para realizar a proxima interação
			pthread_join(b_next_t, NULL);
			pthread_join(t_next_t, NULL);
			pthread_join(p_next_t, NULL);

		}

	}	

		// FINALIZANDO O CALCULO
		// Aqui são feitos os calculos finais para chegar no valor de PI
		
		// pi = a_next + b_next
		mpf_add(pi, a_next, b_next);
		// pi = pi²
		mpf_pow_ui(pi, pi, 2);
		// t_next = t_next * 4
		mpf_mul_ui(t_next, t_next, 4);
		// pi = pi / t_next
		mpf_div(pi, pi, t_next);

		// define o numero de casas de precisão
		int n = 50;
		// imprime PI
		gmp_printf("\n\n   PI = %.*Ff\n\n", n, pi);

	return 0;
}

// Função para calculo de 'a_next' para thread
void *calc_a(void *t){
	// a_next = (a_prev + b_prev) / 2;
	mpf_add(a_next, a_prev, b_prev);
	mpf_div_ui(a_next, a_next, 2);

	pthread_exit(0);
}

// Função para calculo de 'b_next' para thread
void *calc_b(void *t){
	// b_next = (a_prev * b_prev) ^ 1/2
	mpf_mul(b_next, a_prev, b_prev);
	mpf_sqrt(b_next, b_next);
	
	pthread_exit(0);
}

// Função para calculo de 't_next' para thread
void *calc_t(void *t){

	//t_next = t_prev - p_prev (a_prev - a_next) ^ 2
	mpf_sub(t_next, a_prev, a_next);
	mpf_pow_ui(t_next, t_next, 2);
	mpf_mul(t_next, t_next, p_prev);
	mpf_sub(t_next, t_prev, t_next);
	
	pthread_exit(0);
}	

// Função para calculo de 'p_next' para thread
void *calc_p(void *t){
	// p_next = p_prev * 2
	mpf_mul_ui(p_next, p_prev, 2);
	
	pthread_exit(0);
}